Reducing the entanglement growth 

in finite-temperature DMRG transport calculations 



C. Karrasch, J. H. Bardarson, and J. E. Moore 

Department of Physics, University of California, Berkeley, California 95720, USA 
Materials Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 
94720, USA 

Abstract. Finite-temperature transport properties of one-dimensional systems can 
be studied using the time dependent density matrix renormalization group via the 
introduction of auxiliary degrees of freedom which purify the thermal statistical 
operator. We demonstrate how the numerical effort of such calculations is reduced 
when the physical time evolution is augmented by an additional time evolution within 
the auxiliary Hilbert space. Specifically, we explore a variety of integrable and non- 
intcgrablc, gapless and gapped models at temperatures ranging from T = oo down 
to T/bandwidth = 0.05 and study both (i) linear response where (heat and charge) 
transport coefficients are determined by the current-current correlation function and 
(ii) non-equilibrium driven by arbitrary large temperature gradients. The modified 
DMRG algorithm removes an 'artificial' build-up of entanglement between the auxiliary 
and physical degrees of freedom. Thus, longer time scales can be reached. 



PACS numbers: 71.27.+a, 05.60.Gg 
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1. Introduction 

A physical system is usually characterized by its response to perturbations. In transport 
setups, one studies charge or energy currents driven by voltage or temperature gradients. 
From the theoretical point of view this is complicated - computing the quantum 
mechanical time evolution of a system in non-equilibrium is one of the most active 
areas of research in condensed matter physics. When the external perturbations are 
small, one can simplify the problem by resorting to the Kubo formalism. E.g., the 
charge conductivity cr(oj) describes the linear response current J induced by a small 
electric field (we will be more specific below): 

a(cu) ~ J e luJt (J(t)J)dt , (1) 

where the dynamical correlation function (J(t)J) is calculated in thermal equilibrium: 

p-h/t 

(A(t)B) = Tr (p T e im Ae- iHt B) , p T = , (2) 

with H being the Hamiltonian of the system and T denoting the temperature. 
Unfortunately, computing transport coefficients such as c(cu) is generally still difficult: 
Even if one knows the exact thermal density matrix px (or the exact ground state), 
extracting correlation functions, which couple all excitations, remains a formidable task. 

The key question posed and addressed in this work is: How can linear-response 
and non-equilibrium transport properties of one-dimensional (Id) systems at finite 
temperature be calculated efficiently using the density matrix renormalization group 
(DMRG) [HE]? DMRG was originally devised [31 H] as tool to accurately determine 
ground states of Id Hamiltonians. The reason for its success became understandable 
when it was formulated using matrix product states (MPS) [5j El tZJ El E]: The area 
law [10] stipulates that the ground states of Id systems governed by local Hamiltonians 
are only entangled locally; this implies that they can be expressed efficiently using a 
MPS with a small bond dimension x (which encodes the amount of entanglement). The 
MPS describing a given system can be determined variationally - this is the very core 
of a ground state DMRG calculation [31 H] . One way to compute correlation functions 
(A(t)B) at T = is to directly simulate the time evolution (for other approaches see 

Rets, pi nu ca m m na eed 

e~ iHt B (ground state) (3) 

using a time-dependent DMRG framework [TTl HH HH1 EQl EH [22l [231 EH E5] . The 

corresponding algorithm can again be formulated elegantly using matrix product states. 
The physical growth of entanglement implies that the bond dimension needed to 
approximate (A(t)B) to a certain accuracy grows with time. This limits the accessible 
time scales. 

Standard DMRG methods allow computing the time evolution of a pure state 
and are thus not directly applicable at T > 0. Various approaches for simulating 
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finite-temperature dynamics using DMRG have been explored within the literature 
[261 [2TI EH [291 [3Q1 EH EH [331 EH ESJ [361 E71 [381 EH SO]. These include probabilistic 
sampling over an appropriately chosen set of pure states [30] , schemes which time-evolve 
operators instead of states (3TJ [32] , transfer- matrix DMRG [331 Ell ESI [36] , and exact 
representations through purification [371 EH EH]. Purification expresses the thermal 
statistical operator p T as a partial trace over a pure state living in an enlarged 
Hilbert space where auxiliary degrees of freedom Q encode the thermal bath |41j : 

p T = Tr Q | * T )(* T | . (4) 

One of the main advantages of this approach is that all the standard methods for time 
evolving quantum states within DMRG are directly applicable. In this work, we employ 
the time evolving block decimation [T71 122] . 

One of the first applications of finite-temperature 'purification DMRG' to dynamical 
problems was the calculation of spin-spin correlation functions of integrable spin-1/2 
Heisenberg chains [38]. While DMRG yields data which are 'numerically exact' (this 
is verified by comparing with analytic results available for 'non-interacting' models 
[331 EH EH]) , the time scales accessible at finite temperatures are considerably smaller 
than at T = 0. This observation, which might be one reason why studies of finite- T 
dynamics [SU E21 ESI EH ESI EH EE! EH1 SQl S21 SHJ SI] are rarer than their T = 
counterparts, can be understood as follows. Assume we want to compute a ground state 
correlation function, i.e. evaluate Eq. (E|). Under the time evolution, the entanglement 
grows locally around the region on which B acted. In contrast, at T > where one 
needs to calculate (see below for details) 

e- im B\^ T ) , (5) 

the entanglement generally increases homogeneously throughout the whole system. This 
holds true even for B = 1, i.e. when the state \^t) is exposed to a supposedly trivial 
time evolution [10]. The reason for this is that purification is not unique and various 
representations of the same px differ in their degree of entanglement between the 
auxiliary and physical degrees of freedom. Even when one starts out with a purification 
that minimizes this entanglement, it rapidly grows under the DMRG time evolution. 
This observation naturally leads to the question: Can that very same 'gauge freedom' 
be used to our advantage to undo the non-physical growth of entanglement? In other 
words, is there a unitary transformation Uq(t) acting on the auxiliary Hilbert space 
Q which removes the artificial entanglement? In Ref. [10] we proposed to time-evolve 
Q backwards in time with the physical Hamiltonian acting on the auxiliary degrees of 
freedom: 

U Q (t) = e +lflt . (6) 

This renders the time evolution of \^t) trivial, and the evaluation of Eq. ([5]) is therefore 
eventually only plagued by an entanglement building up around the region where B acts 
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in complete analogy to the ground state calculation (the physical reason being quasi- 
locality [311 [32]). This generically leads to a slower increase of the bond dimension x, 
and thus longer time scales can be reached. 

In Ref. (40], the potential of the modified DMRG algorithm was demonstrated for 
the spin-spin correlation function (S^^S^) of the XXZ chain at zero anisotropy A = 0, 
which maps to free fermions and thus allows for an exact (benchmark) solution. A more 
thorough comparison between the numerical effort of the standard [38] and modified 
[40] algorithms in calculating {S^{t)S^) for arbitrary A can be found in Ref. [32], where 
further optimizations using operator-space DMRG were explored. One of the ideas of 
Ref. [32] is to use time translation invariance to rewrite Eq. (j2J) as 

(A(t)B) = (A(t/2)B(-t/2)) , (7) 

which allows accessing times twice as large without additional effort. 

From the point of view of physical applications, the modified DMRG schemes were 
used to investigate current correlation functions and the Drude weight of the integrable 
XXZ chain at intermediate to high temperatures [401 E], scaling properties of the DC 
conductivity in presence of non-integrable perturbations [12], non-equilibrium induced 
by temperature gradients [43], and spectral functions of hardcore bosons [31]. For 
reasons of completeness, we mention other approaches to linear response transport 
properties of the XXZ model (and related ones). These include exact Bethe ansatz 
calculations [471 SHI SHI ED], integrability arguments [SH EZl E31 E3 ESI EE], field 
theories [351 ESI EH EH EH], quantum Monte Carlo [601 EH E2], exact diagonalization 
[631 E3 ESI ESI EZ], transfer matrix DMRG [351 EE], and dynamical DMRG [H]. 
For a non-exhaustive list of prior works on non-equilibrium thermal transport see 
Refs. [23lESEHl[TQl[7J[7j[7iil[^ 

The first and foremost goal of this paper is to present extensive quantitative data 
on how the build-up of entanglement is reduced if the auxiliary Hilbert space Q is 
time-evolved with the physical Hamiltonian but reversed time. We focus on transport 
properties in linear response [40 1 142 1 144] and in thermal non-equilibrium [43]. Specifically, 
we study (i) homogeneous gap less and gapped spin- 1/2 XXZ chains, also in presence of 
various perturbations which break integrability, (ii) the quantum Ising model, and (iii) 
impurity setups of a quantum dot connected to non-interacting leads. For temperatures 
from T = oo down to T/bandwidth = 0.05, the modified DMRG algorithm leads to a 
slower increase of the MPS dimension \- Only at low T does the standard approach 
Uq = 1 become more efficient (see Ref. [32] for details). For most (but not all; we 
will be more specific below) problems studied in this work, T/bandwidth = 0.05 is a 
low enough temperature to correspond to the T = limit, which one can establish, 
e.g., by comparing with field theory results [77]. We show the typical behavior of x on 
the relevant physical time scales for each problem at hand (e.g., on the time scale at 
which non-equilibrium currents generically reach their steady state). We reiterate how 
Eq. (JTJ) can be implemented within the purification approach and illustrate (following 
Refs. (SUES]) how it allows accessing times twice as large. 
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As the second purpose of this paper, we present technical details of the 
implementation of the algorithm. E.g., we show how to time-evolve next-nearest 
neighbors (closely following Ref. [2]), which is necessary because physical degrees of 
freedom are separated by an auxiliary site due to the purification. We discuss the 
numerical accuracy of our data, compare with exact results, investigate to what extend 
the choice of Uq{t) = exp(iHt) is optimal [SUES], and provide further evidence for the 
reliability of linear prediction extrapolation schemes [38| W2\ H5| 06] for the spin-spin 
correlation function of the isotropic Heisenberg chain. 



2. Models and methods 



2.1. Models 

As a prototypical model, we consider a chain of interacting spin- 1/2 degrees of freedom 
gx,y,z governed by local Hamiltonians 

or equivalently spinless fermions through a Jordan- Wigner transformation. By choosing 
the couplings J n , A n , and b n appropriately: 

fl n odd , (-l) n b 

J n = { , , A n = A , b n = , (9) 

I A n even * 

we can study systems which are gapless or gapped, and investigate the role of 
integrability. For A = 1 and 6 = 0, Eq. (jHJ) can be diagonalized via Bethe ansatz 
|79j ; the model is non-integrable otherwise. The energy spectrum is gapless for |A| < 1 
and gapped for A > 1. A gap opens for A < A c or b > b c , where A c < 1 and b c > if 
— 1 < A < —l/y/2 [HU EH 132] . In addition, we study the quantum Ising model 

h n = -AS z n S z n+l -g{S* + Sl +l ) . (10) 



2.2. Transport properties 

2.2.1. Linear response. We first consider a homogenous system of size L governed 
by H — X^n=i within linear response. The optical charge (C) and energy (E) 
conductivities can be computed via the Kubo formula 

i roc i _ --w/T poo 

°{") = -J / e^{[J(t), J])dt = — / e^(J(t)J)dt , (11) 

where the corresponding current operators J = j n are defined through a continuity 
equation [82"] : 

L-l 

dth n = jE,n. — jE,n+l =>• <7e = i / ][h n -l, h n ] , 

Zl ( 12 ) 



dtSn — JC,n — JC,n+l => Jc — i ^[^n-l> Sft 



n=2 
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Figure 1. (Color online) The non-equilibrium setups studied in this work. A: 
Two interacting chains of length L/2 which are initially in thermal equilibrium at 
temperatures Tl and Tr are coupled at time t = 0. B: Two non- interacting chains are 
coupled at time t = via an interacting resonant level model (IRLM). 



One usually decomposes the real part of as 

Re ct(uj) = 2nD5{oj) + <T reg (c<;) , (13) 

where the so-called Drude weight D is the prefactor of the singular contribution. D can 
be determined from the long-time asymptote of the current-current correlation function 

D = hm hm \ K J ' . (14) 

t— >oo L— >oo 2LT 

As already stated above, a time- dependent equilibrium correlation function is defined 

as 

P -H/T 

(A(t)B) = Tr (p T e iHt Ae- iHt B) , p T = —- , (15) 
with Z T = Trexp(— H/T) denoting the partition function. 

2.2.2. Non- equilibrium. In addition to linear response, we study two thermal non- 
equilibrium setups. The first [labeled 'non-equilibrium A' and depicted in Figure ED(A)] 
is introduced via the following protocol: We initially consider two separate chains, 

1/2-1 L-l 

H = H L + H R = h n+ Yl h "> ( 16 ) 

n=l n=L/2+l 

each being in thermal (grand-canonical) equilibrium at temperatures Tl and Tr. The 
corresponding density matrix factorizes, 

exp(-Hi/Ti) -rr, n>7\ 
Pq = Pl®Pr, Pi = Trexp (_ffy T .) » i = L,R. (17) 

At time t — 0, the chains are coupled through fiLii, and the time evolution of any 
observable A is computed using H = H + h L / 2 : 

(A(t)) = Tr p(t)A , p{t) = e iHt p e- iHt . (18) 

In the second setup ['non-equilibrium B'; see Figure HB)] we investigate two non- 
interacting chains A = b = 0, A = lof length L/2, 

L/2-1 L 

H = H L + H R = h n + J2 h n> ( 19 ) 

n=l n=L/2+2 



Reducing the entanglement growth in finite-temperature DMRG 7 

at different temperatures Tl and Tr. At t = 0, they are coupled via an interacting 
resonant level model: 

friRLM = t'{S x L/2 S x L / 2+1 + S y L/2 S y L/2+l + US z L/2 S z L/2+l ) 

+ ' {^1/2+1^1/2+2 + ^L/2+1^1/2+2 + ^ ^1/2+1^1/2+2) • 

The site n = L/2 + 1 is initially in an equal superposition of up and down states (i.e., 
formally at infinite temperature). 

2.3. Purification. 

One way to evaluate Eqs. f JT5|) and ffT8j) within the density matrix renormalization group 
is to purify [JT] the thermal density matrix by introducing an auxiliary Hilbert space 
Q: 

e -H/T 

p T = ^^ = Tr Q \^ T ){^ T \ . (21) 

At infinite temperature, where px is proportional to the unit operator, one can 
analytically write down the purification: 



P 



1 1 L 

^ = ^=Tt Q \^o )^o \=Tt q H\^J(^J . (22) 



n=l 



In order to exploit conservation laws within the DMRG algorithm, it is convenient to 
choose 

|v&« ) = — (\a n =ta nQ =i) - \a n =i,a nQ =t)) , (23) 
where <r n =t, 4- denotes the eigenbasis of S*, and Q is spanned by 

Q = span{|(T lQ ...(T LQ )} . (24) 
One readily verifies that indeed 

The finite-temperature purified state |\l/x) is obtained from l^oo) by applying an 
imaginary time evolution [2], 

Pt = ^ = ^Trge-^|*oo)(*oo|e-^ 2T 

Zt Zt (26) 

where we formally replaced 

H^H®l Q . (27) 
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The correlation function of Eq. (j!5p is exactly recast as 

(A(t)B) = g ^ , (28) 

(^oole^-^oo) 

and for two time arguments one similarly finds 

(yDrpl piHtA -iHt iHt' d -iHt' I q, \ 

= (Wk) ' (29) 

In order to evaluate Eq. (J7|) for A = A^" = i?, it is sufficient to compute 

The expectation value of any observable A with respect to the time-dependent (non- 
equilibrium) density matrix of Eq. ( !T8|) is given by 

,a, A \ rp ^ (y 00 \e-&-^e im Ae- iHt e-&-*%\V 00 ) 

(A(t)) = Tr P (t)A = • ( 31 ) 

24. DMRG algorithm 

2. 4-1- Initial state. It is convenient to first express the initial state |^oo) in terms of a 
matrix product state [5j [6j [7J |8] , 

= ^ A^A^Q . . . A° L A aL Q |a lf 7 lQ . . . a L a LQ ) , (32) 

where 

^aia i+ i ~~ ^-aj- \ia i+1 • (33) 

Here and in the following <ji is a short hand for either a physical or auxiliary degrees of 
freedom: Ui G {c n ,a nQ }. The initial matrices corresponding to Eq. ( 123]) read 

r CT - =t = (i o) r CTn =i = (o - 1) 

r CT "« =t = (o i/\/2) T r CT -Q =J - = (i/y/2 o) T , (34) 

as well as A 1 = 1. Normalization generally stipulates 

V(A*„C ) t (A!,r <Jl , ) = 6 Q .w (35) 

as well as 

E ( r ^ +1 A 2i)( r <a l+1 A 2 + 1 1 ) t = • ( 36 ) 
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2.4-2. Trotter decomposition. In order to successively apply the imaginary and real 
time evolutions operators 

e~ XH = e~ AXH e- AXH ... , A = AA + AA + . . . (37) 

to the initial state l^oc), we factorize them either via a second or a fourth order Trotter 
decomposition [2]. The second order scheme reads 

e -AXH ^ e -AX/2H le -AXH 2e -AX/2 Hl ? (3g) 

where the H = Hi + H 2 , and H 12 contain only terms which mutually commute: 

H X =Y,K, H 2 = J2 h n- (39) 

n odd n even 

In order to reduce the error, one can employ a fourth order decomposition [M] , 

e~ AXH « f/(AAx) 2 U (AA 2 ) £/(AAi) 2 , (40) 

where 

C/(AAj) = e -AA i ff 1 /2 e -AA i if 2e -AA i fr 1 /2 ? (41) 

and 

AAi = 4 _ 41/3 AA , AA 2 = AA - 4AAx . (42) 

We typically employ a second order Trotter decomposition with a time step of A/3 = 
0.005 during the imaginary time evolution from = down to the physical temperature 
(3 = 1/T and a fourth order decomposition with a time step of At = 0.2 during the 
real time evolution. This is sufficient for all problems studied in this paper (which we 
verify by comparing against data obtained for smaller A/3 = 0.0025 and At = 0.1). The 
Trotter decomposition is typically not a significant source of error within time-dependent 
DMRG calculations |68|. 



2.4.3. DMRG step. The physical Hamiltonians of Eqs. (IE]) and ( TTOl) couple matrices 
with indices a n and cr n +i which in the matrix product state of Eq. (|32|) are next-nearest 
neighbors - they are separated by a matrix associated with an auxiliary degree of 
freedom a nq . Through Eq. (T27j), the local Hamiltonians of Eqs. ([8]) and f lTUj) are formally 
replaced by 

K -> h n ^ Wn Q )(o- nQ \ ■ (43) 

Thus, after each application of exp(— A\h n ), two singular value decompositions are 
carried out to update three consecutive matrices; to keep the notation transparent, let 
us denote those matrices as: 

p CTl y jVi A i y jl 

T «2 y f<T2 A 2 y ^2 ^ 
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and furthermore abbreviate h = h n . This 'DMRG update' can be achieved by a simple 
and straightforward generalization of the protocol for nearest neighbors outlined in 
Ref. [2]. First, one forms the three-site tensor, 

^<T a lC T 2 _ AO p(70 A 1 pO-1 A 2 p0 2 A 3 (Ar\ 

1Q13 «0 ffio«l «1 «1<12 12 «2«3 "3 ' V^"/ 

aia2 

which is then acted on by exp(— AXh): 

1013 / / 1013 lo OiO 2 ' V ' 

CT 1°"2 CT 3 

Now, two singular value decompositions (SVDs) are applied to appropriately reshaped 
tensors in order to obtain the new matrices A* and 



0013 ^ (10 oq J, (13 0102 J 



/]^(aogtt),ttiAft 1 (V ,t ) ttl , 



(0-10-213) 



/ , 10 ^ no < 

11 - „ v— 

f CT 
1 a a l 



(47) 



00 

1011 ' 



and likewise 



^11,(010-213) ^(1101), (0203) 



12,(0213) 



02 



VA 1 A 1 - 1 !/"" 1 A 2 (V+r 2 a^a 3 

/ j 11 11 1112 12 \ /1213 13 a 



(4? 



L 3 

13 

12 



p°l p°2 
la l«2 L a 2 a 3 



In case of a real time evolution, the normalization conditions of Eqs. (155|) and 
are preserved automatically (up to errors associated with a potential truncation of the 
matrices; see below). For a non-unitary exp(— A/3h), however, the matrix product state 
needs to be brought back to its canonical form in order to simplify the evaluation of 
expectation values, and, more importantly, because otherwise the truncation below is 
not optimal. Normalization can be achieved approximately by successively acting with 
unit operators AA = (22] or exactly using the procedure outlined in Ref. [83] . 



2.4-4- The discarded weight. Up to this point, the presented DMRG algorithm to 
evaluate Eqs. (j2"8|) -( lHT|) is exact (except for the Trotter error which can be made 
sufficiently small for all problems studied in this work). However, the dimension of 
the matrices in Eq. ( )32l) increases during each of the singular value decompositions in 
Eqs. ( flT)) and (HSj) . generically by a factor of two (corresponding to the number of local 
degrees of freedom). For brevity, let us neglect that for open boundary conditions the 
matrix dimension is position-dependent and simply denote it by \ (a more thorough 
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discussion can be found in Ref. [2]). The computational effort of the DMRG algorithm 
is dominated by the cost of the singular value decompositions, which scales as x 3 - 
Increasing \ — > 2\ during each sub-step of a Trotter step is numerically unfeasible (and 
fortunately generically unnecessary). Thus, we truncate the matrices to a dimension 
X' < %X by dropping the indices a* which correspond to the smallest singular values A ai . 
For a normalized state [i.e., if Eqs. fl35|) and f )36|) hold], this is the best approximation 
with respect to the Hilbert space norm. The error (i.e., the norm-distance of the states 
before and after truncation) is given by the sum of all squared singular values which are 
discarded during all sub-steps of one Trotter step. This so-called discarded weight is 
the key numerical control parameter; the DMRG algorithm becomes exact as it is sent 
to zero. 

One can strictly enforce a fixed discarded weight by analyzing the singular values 
during all SVDs carried out within one Trotter substep and then retrospectively choosing 
X 1 appropriately (or by simply repeating the substep if x' > x)- More pragmatically, 
one can fix x' — X an d increase x by A% after one substep if the total discarded 
weight exceeds a pre-defined threshold e. While this simple approach certainly becomes 
inapplicable if one aims at strictly maintaining a discarded weight that is very small, it 
works for all problems studied in this paper, with Ax ~ 10—20 being a reasonable choice. 
We observe that the actual discarded weight is smaller than 2e for all data presented in 
this paper [it is imperative to always monitor the real discarded weight which determines 
the actual error; this is illustrated in the inset to Figure |3]^a)]. In the following we neglect 
this difference and simply refer to our control parameter e as 'discarded weight'. We 
successively lower e from 10~ 8 — 10~ 13 during the imaginary time evolution and from 
10 -3 — 10 -7 during the real time evolution, which turns out to be sufficient to obtain 
'accurate' results for all physical quantities of interest. 'Accurate' means 'sufficiently 
converged with respect to lowering e', or 'in reasonable agreement with exact data' in 
case that the problem allows for an analytic solution; all this is illustrated in Figures 
|3]^a), |5]^a), |6]^a), and|8]^a) and will be discussed in detail below. We stop our simulations 
once the matrix dimension has reached values around x ~ 1500 — 2500. 

2.5. Evolve the auxiliaries! 

An exact modification to the finite-temperature DMRG algorithm (which allows 
accessing longer time scales) was recently introduced in Ref. [ID]: One has the analytic 
freedom to apply any time- dependent unitary transformation 



to the in principle inert auxiliary sites a nQ - physical quantities are determined by the 
trace over Q and are thus not affected by UQ(t): 




Q 



Q 



Q 



(49) 



[U Q (t),* n ] = 0. 



(50) 
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For example, Eq. (130 j) can be rewritten as 

(y T \e im Ae- im e im 'Be lHt '\y T ) = (m T \U\t)AU(t)U\t')BU{t')\^ T ) , (51) 

where 

U(t) = e- mt U Q (t) . (52) 
Put differently, purification is not unique. It turned out [311 [32j 001 021 03] that choosing 

{7q(£) = e +lHt , if = H(a n — > a nQ , magnetic fields reversed) , (53) 

i.e., time-evolving the auxiliary sites with the physical Hamiltonian but reversed time, 
leads to a slower increase of \ and thus longer time scales can be reached. An intuitive 
way of understanding this will be given below. Only at low temperatures does using 
Uq = 1 become more efficient [3H[32]. In practice it is most convenient to implement 

-iHtJHt _ -iH At AH At -iH At AH At (r A \ 



3. Results 



In this section we present our main results. We first quantify how the entanglement 
growth is reduced by time-evolving the auxiliaries with the physical Hamiltonian but 
reversed time. In an extension of earlier studies (see Refs. [311 [321 HO] and the discussion 
in Sec. [Q, we study gapless and gapped XXZ chains in presence of non-integrable 
perturbations and focus on transport problems both in linear response and out of 
equilibrium. We investigate to what extend £/ aux (t) = exp(itH) is optimal within our 
approach. In passing, we provide another example of the stability of linear prediction 
(see also Refs. [3 M, M, S3 06]) and recast (one of) the ideas of Refs. [SH [32] in the 
language of purification. 

3.1. Prelude: Time evolution of \^t) 

It is instructive to first consider the trivial case of A = B = 1 in Eq. ([28]) . i.e. 

x = (* r |e""e-""|g r ) = ^T\e' H %{t)e^ m U Q {t)\m T ) 

Figure [2] shows the evolution of the bond dimension x during the calculation of 
e~ Uq{C)\&t)- For the standard approach Uq{t) — 1, % increases with time: the 
state e~ lHt \^T) which purifies the density matrix becomes successively more entangled. 
Choosing Uq{£) = e lHt removes this artifact. This immediately indicates why longer 
time scales can be reached in DMRG evaluations of, e.g., current correlation functions 
e~ UQ(t)j n \i&T) ■ The entanglement only builds up locally around site n (the physical 
reason being quasi- locality [3T1 132]). and likewise for our non-equilibrium setup it grows 
locally around the interface region over which the initial temperature gradient falls off. 
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Figure 2. (Color online) MPS dimension during the time evolution of the state |^t) 
which purifies the thermal density matrix of a XXZ chain. 



The artificial global build-up of entanglement is removed if the auxiliaries are evolved 
in time. Only at low T choosing Uq = 1 becomes more efficient [HTJ 132] since the time 
evolution of the ground state is trivial (the latter is indicated in Figure [2j x grows more 
slowly at low T = 0.2 than it does at high T). 

An intuitive understanding of why the particular choice of Uq(t) = e +lHt removes 
the 'artifical' entanglement growth was recently provided in Refs. [31j [32]: In an 
operator-space formulation, the modified DMRG algorithm corresponds to a Heisenberg 
time evolution of the matrix product operators representing A and B in Eq. (I5T1) . If A 
(and likewise for B) is local, most terms in the first Trotter step e *# A *^4 e -^ A * cancel 
out. As time evolves, 

...e iHAt e iHAt Ae- iHAt e- iHAt ... , (56) 

entanglement can only build up gradually around the region on which A acted due 
to the so-called light-cone effect in nonrelativistic systems: Lieb-Robinson bounds [85] 
generically stipulate that correlation functions (A(t)B) of local operators A and B fall 
off exponentially for x > vt, with x denoting the spatial distance between the regions 
on which A and B act, and v being some velocity with which excitations propagate. 



3.2. Reduction of the growth of entanglement: N on- equilibrium 

In this section we illustrate quantitatively the effects of time-evolving the auxiliaries for 
the non-equilibrium setups sketched in Figure [U We start by studying an interacting 
XXZ chain (A ^ 0) with two additional perturbations (dimerization A < 1 and a 
staggered field b > 0) that both render the system non-integrable [see Eq. ([8]) for a 
precise definition]. At time t = 0, two 'semi-infinite' chains (typical lengths being 
L/2 ~ 100 — 200) prepared in thermal equilibrium at temperatures T^jj are coupled 
by hi/2 to an overall 'translationally-invariant' chain. The temperature gradient drives 
an energy current (jE,n(^)) whose typical behavior is shown in Figure [3]^a). For 6 = 
the current at the interface n = L/2 saturates to a finite value on a scale t ~ 1 — 10 
irrespective of whether the system is gapless or gapped (indicating that either the non- 
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t 



Figure 3. (Color online) DMRG time evolution of an XXZ chain of length L = 200 
featuring an initial sharp temperature gradient Tjr, ^= Tr. The system is gapless for z- 
anisotropics |A| < 1 and gapped otherwise; it becomes non-integrable in presence of a 
finite dimcrization A < 1. Note that for A = 0.5 and A = 1 asymptotic low-temperature 
behavior described by a field theory [77] sets in for T < 0.2. (a) Energy current 
at the interface. Choosing a discarded weight control parameter around e ~ 10~ 6 
at a Trotter stepsize of At = 0.2 is sufficient to accurately obtain its steady-state 
value. Inset: The actual discarded weights during a whole Trotter step and during 
one substep for A\ = 20. The difference to e is explained in the main text, (b) 
Evolution of the dimension of the corresponding matrix product state. If the auxiliary 
degrees of freedom which purify the thermal density matrix are time-evolved with 
the physical Hamiltonian but reversed time (which is an exact modification to the 
standard algorithm) , the build-up of entanglement is reduced. The computational cost 
of a singular value decomposition (which dominates the DMRG algorithm) scales as 
X 3 - 



integrable dimerized chain supports dissipationless transport or that its current decays 
on a hidden large scale; see Ref. [13] for further details). The steady-state current of 
the XXZ chain is of the simple 'black-body' form [H2 [76J [77J, EE] 

}im{j E>n (t)) = f(T L )-f(T R ) , (57) 

implying that it is determined solely by the linear conductance G(T) ~ drf{T) for 
arbitrary large Tl — Tr. 

lim0E,n(*))~ / R G{T)alT . (58) 

t^oo J Tl 

This agrees with a field theory result [77] at low temperatures; asymptotic low-T 
behavior sets in around T < 0.2. 

The time evolution of the corresponding bond dimension \ is shown in Figure 
[3]^b) for the integrable XXZ chain both for parameters where it is gapless (A = 0.5) and 
gapped (A = 2). Figure|H(a) shows the same in presence of non-integrable perturbations 
as well as for the quantum Ising model of Eq. ( flOl) at criticality g = 1. In all cases and for 
all temperatures from T^r = oo down to Tl_r = 0.1 (which for this problem corresponds 
to zero temperature), time evolution of the auxiliaries leads to a slower increase of x; 
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Figure 4. (Color online) (a) The same as in Figure [3jb) but in presence of two 
perturbations (dimerization A < 1 and a staggered magnetic field b > 0) rendering the 
model non-integrable. Inset: Quantum Ising model at criticality 3 = 1. (b) Dimension 
of the MPS during the time evolution of two non-interacting XXZ chains of lengths 
L/2 = 100 which feature different temperatures Tl and Tr and are coupled at time 
t = via an interacting resonant level model. 

note that the computational cost of a singular value decomposition (which dominates 
the whole DMRG algorithm) scales as x 3 - 

The same reduction of x holds for the non-equilibrium impurity problem depicted 
in Figure H(b) where an interacting resonant level model defined in Eq. (1201) is coupled to 
two non- interacting leads (A = b = 0, A = 1) at different temperatures. An exemplary 
evolution of the MPS dimension is shown in Figure H(b); a discussion of the (transport) 
physics at small interactions or zero temperature can be found in Refs. (87]. Note that 
due to the appearance of a Kondo-like scale, T/bandwidth = 0.05 does not necessarily 
correspond to zero temperature for this problem |87j . 

The discarded weight e = 10~ 6 in Figure [3]^b) is chosen small enough to 'accurately' 
determine the current and in particular its steady-state value. This is illustrated in 
Figure |3^a) where e is successively lowered from e = 10~ 4 to e = 10~ 7 . Moreover, the 
non-interacting case A = b = 0, A = 1 allows for an exact evaluation of the steady- 
state current [76], which can be used to benchmark the accuracy of the DMRG data (a 
comparison can be found in Ref. [12]). The inset to Figure [3]^a) shows a generic example 
for the actual discarded weight (for the precise definition of e see Sec. 12.4.41) . It is always 
bound by 2e; the total discarded weight during a complete Trotter step e - lHAt e iHAt j s 
typically bound by lOe. 

The discarded weight determines the norm distance of the exact time-evolved state 
and its DMRG approximation; this holds for any choice of Ug(t). However, for Ug(t) = 1 
(i.e., the standard approach of inert auxiliaries) entanglement builds up artificially even 
far away from the interface n = L/2 where the initial temperature gradient is applied, 
and the corresponding truncation error cannot immediately manifest in local observables 
(such as the currents around n ~ L/2). This is again due to the Lieb- Robinson bounds 
discussed in Sec. 12.51 loosely speaking, in our case A and B can be thought of as the 
current jLji an d as a local time evolution operator e~ lAthn ^ L / 2 , respectively. It is thus 



Reducing the entanglement growth in finite-temperature DMRG 



16 



A 
O 
—} 

~~0 
—} 

V 

CD 



0.1- 

0.05- 

0- 

0.1- 

0.05- 



T=0.5 




(a) 



X-70 
X-280 

JZioSn fixed x=200 fixed x=200 

* 1JUU (standard) (evolve aux) 

\... 



1/ 



e=1Cr 4 

"'~e=i(r ! 



X~90 
X-330 

X~750 
X-1500 

e=10" 



T=2 



E=10" 4 



10 



20 




Figure 5. (Color online) (a) DMRG calculation of the global charge current 
correlation function ( Jc(t)Jc) of the integrable XXZ chain (length L = 100) in thermal 
equilibrium. The calculation is stopped once x reaches values around x ~ 1500; the 
numbers in the plot denote x & t time t = 12 for the different discarded weights (the 
auxiliaries arc time-evolved), (b) Evolution of the corresponding MPS dimension for 
different anisotropics A = 0.71 (gapless) and A = 3 (gapped). 



instructive to carry out a calculation at a fixed MPS dimension x rather than a fixed 
discarded weight. The result is also shown in Figure EJ^a) - it clearly rules out the 
possibility of the artificial error 'propagating' so slowly that it would not manifest on 
the relevant physical time scale t ~ 10. 

3.3. Reduction of the growth of entanglement: Equilibrium 

We now turn to current correlation functions in thermal equilibrium. Following Eqs. (|T2|) 
and ( )28l) . we need to evaluate 

e- lHt U Q (t) 3n \y T ) . (59) 

The AC conductivity c(uj) is determined by the Fourier transform of (J(t)J) via 
Eq. ( TTTj) . The integrable gapless XXZ chain, on which we focus in this section, supports 
dissipationless transport at finite temperature, i.e. its Drude weight D in Eq. ( IT3|) is 
finite. D can be obtained from the long-time limit of the current- current correlation 
function via Eq. ( fl4l) . The relevant time scale in this problem is thus the scale on which 
(J(t)J) saturate to their asymptote. 

DMRG data for the charge current correlation function at A = 0.71 and for 
various discarded weights is shown in Figure |5]^a). At intermediate to high temperatures 
T > 0.5, we can access time scales on which (J(t)J) saturates [we stop our simulation 
once the MPS dimension has reached values around x ~ 1500; see the numbers given 
in Figure GJa)]. More details can be found in Refs. [101 03]. A fairly large e ~ 10~ 5 is 
sufficient to obtain D quantitatively. This is supported by the comparison with the exact 
solution for A = shown in the inset to Figure|6^a). Note that (Jc(t)Jc) is constant for 
A = since Jc is conserved; however, every single term (jc,n(t)jc,m) that contributes 
to the sum is time-dependent. It is again instructive to carry out a calculation using 
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Figure 6. (Color online) (a) Equilibrium spin and global charge current correlation 
functions S zz (k,t) = J2 n ^^n+L/zi^ln) and (J C (t)J c ) of the XXZ chain at 
A = 0. DMRG results are compared with the exact solution obtained by mapping the 
model to free fermions. Since Jq is conserved, (Jc{t)Jc) is constant (we dropped some 
DMRG data points to improve visibility). However, every single term (jc,n(t)jc,m) 
that contributes to the sum is time-dependent, (b) Demonstration of the stability of 
linear prediction (see also Ref. [38] and Refs. [2] [31] [42] [45] [46] ) for the isotropic chain 
A = 1. Dashed lines were obtained by fitting the DMRG data for times tat G [2-5, 5] 
and subsequent extrapolation to t > 5 using linear prediction. The extrapolated data 
almost coincides with the DMRG data for t > 5 (solid lines). 



a constant x rather than a fixed discarded weight; the result is shown in the upper 
panel of Figure Efa), illustrating that for the standard approach Uq{€) = 1 the 'artificial 
errors' propagate so fast that they certainly manifest on the relevant time scale of this 
problem. 

For all temperatures from T = oo down to T = 0.125 time evolution of the 
auxiliaries leads to a slower increase of the dimension of the matrix product state. 
This is depicted in Figure Efb). One might expect low-T behavior of D(T) to set in for 
T = 0.125 (see Refs. [351 1361 E]), but the current correlation function decays so slowly 
that we cannot access its asymptotics without extrapolation. Exemplary data for the 
gapped phase is shown in the lower inset to Figure [5^b). 

3.4. Spin correlation function 

We briefly present data for the spin-spin correlation function 

S"(n, t) = (S z n+L/2 (t)S z L/2 ) , S zz (k, t) = Y, e ikn S**(n, t) , (60) 

n 

which provides a standard testing ground for DMRG approaches to time-dependent 
correlation functions at finite temperature [33]. 13^1 138] HO] . The XX chain A = again 
allows for an exact solution by mapping the model to free fermions. Figure Efa) shows 
a comparison of DMRG data with the exact result for T = 1 and T = 0.1 (which 
for this situation corresponds to low temperature). A large discarded weight e ~ 10 -5 
is sufficient to reproduce the analytic result up to times t ~ 30. If the auxiliaries 
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11 

Figure 7. (Color online) Entanglement entropy of the state \^>t) after time-evolving 
it up to a time t = 4. The physical degrees of freedom are time-evolved via exp(— iHt); 
the auxiliaries are time-evolved using Uq(t) = cxp(+iHt + itr/N) where N is given by 
Eq. (|63|) (solid lines) and Eq. (|64j) (dashed lines). 

are time-evolved, the MPS dimension increases only moderately to x ~ 100 — 200 in 
contrast to the standard approach (see Ref. [10] and also Refs. [33| IM] for transfer- 
matrix DMRG data). For finite A, the choice of Uq{t) = e lHt still outperforms Uq = 1 
(expect at low temperatures), but the difference becomes successively less pronounced 
[32] . For the isotropic chain A = 1 and all T > 0.1, the accessible time scales are 
about a factor 1.5 — 2 larger if the auxiliaries are evolved in time. In Ref. [38] (see 
also Refs. [21 [311 H21 H3 HE]), linear prediction was first established as an accurate 
tool to extrapolate correlation functions by considering exactly-solvable models and 
subsequently used to compute the Fourier transform of S zz (k,t) of the isotropic chain. 
For reasons of completeness, we revisit the isotropic chain and compare the result of 
linear prediction with DMRG data for larger times. Figure Mjo) illustrates that the 
agreement is good. 

3. 5. Optimization of Uq 

Finally, we shortly investigate to what extend our choice of e lHt for the time evolution 
of the auxiliaries is optimal (more details on the optimization of Uq can be found in 
Refs. [3H[32]). To this end, we consider 

U*(t) = e +iSt+itr > N , (61) 

with N begin an arbitrary Hermitian matrix. We compute the entanglement entropy 

S cnt (t, V) = "2 B A - /2 ) 2 A t /2 " 2 B A ^ e ) 2 lo S 2 A «ll Q (62) 

between the left and right halves of the time-evolved state e~ Uq(t)\^T) which purifies 
the density matrix. We have verified that for arbitrary N coupling nearest neighbors, 
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Figure 8. (Color online) (a) Equilibrium spin correlation functions S zz (n,t) = 
°f the isotropic XXZ chain for various discarded weights. The 
calculations were stopped once the MPS dimension reached values around \ ~ 2000. 
(b) By rewriting S zz {n = 0,i) = (S z L/2 (t/2)S z L/2 (-t/2)) and time-evolving S z L/2 (±t/2) 
independently the time scale t can be accessed by a MPS of a smaller dimension \ 
(this idea was introduced in Refs. p H 132]). 



Sent{t,v) is a ^ least quadratic in 77. Figure [7] illustrates this for the two choices 

L-l 

* = 1 + 2 E + ^ Q+1 ) (63) 

as well as 

L-l 

^ = E = ^ (® n Q =^ a n Q +i =t I • (64) 

nQ=l 

As mentioned above, Refs. [3T| [32] present thorough details on how to further optimize 
finite-temperature dynamical DMRG. We here reiterate one of the main ideas using the 
language of purification. In thermal equilibrium, one can recast any correlation function 
exploiting time translation invariance: 

(A(t)B) = (A(t/2)B(-t/2)) . (65) 

Thus, one only needs to carry out time evolutions up to times \t/2\ in order to compute 
(A(t)B). If A = = B (e.g., for current-current correlation functions), it is sufficient 
to perform a single calculation 

e lHtl2 U ] Q (t/2)Ae- lHtl2 U Q {t/2)\m T ) . (66) 

This implies that one can access a time scale twice as large without much additional 
effort. This is illustrated in Figure HD[b) for the spin- autocorrelation function of the 
isotropic Heisenberg chain. 



4. Summary 

We presented extensive quantitative data for how the growth of entanglement in 
finite-temperature dynamical DMRG calculations can be reduced by time-evolving the 
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auxiliary degrees of freedom which purify the thermal statistical operator. The time 
evolution of the auxiliaries is an exact modification to the standard algorithm. Our 
particular focus was on energy and charge transport properties both in linear response 
(described by current-current correlation functions) and in non-equilibrium driven by 
temperature gradients. We studied a variety of gapless and gapped integrable and 
non-integrable spin- 1/2 chains (i.e., interacting spinless fermions), the quantum Ising 
model at criticality, and impurity (quantum dot) setups. For all temperatures from 
T = oo down to T/bandwidth = 0.05, which for most problems investigated in this 
work corresponds to low temperatures, the build-up of entanglement is reduced. This 
speeds up numerics and allows to eventually access longer time scales. 
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